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I. Introduction 

A three-dimensional, finite volume, time accurate, TVD flow code has been developed 
recently to solve the f ull compressible Navier-Stokes equations. 1 This flow solver employs 
Roe’s flux difference splitting for the inviscid fluxes. A unified formulation was introduced 
for the interpolation of the fluxes on the cell surfaces so that central or upwind schemes of 
various orders of accuracy can be selected by changing parameters in the formula. A second 
order predictor corrector scheme is used for the time integration of the equations. In Ref. 1 
we have shown that the scheme has very good shock capturing capability. 

A moving shock impinging on a circular cylinder produces interesting shock deflection 
phenomena. As the shock move along the solid surface, shock reflection changes from sim- 
ple reflection to Mach reflection, causing complicated shock and slip line patterns. These 
phenomena have been studied both experimentally and numerically. 2,3 

In the present work, the above described computer code is used to study the problem 
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of moving shocks passing objects. This problem not only requires the scheme to have good 
shock capturing capability (non-oscillatory) and good time accuracy, but also requires the 
scheme to be able to treat regions of very low Mach number, for the flowfield in front of the 
shock has zero Mach number and near the stagnation point in front of the object f e g. a 
circular cylinder) the Mach number is very small. Codes that use ADI schemes are known 
to have difficulties in treating such regions. 

The purposes of this paper is to demonstrate the ability of the present numerical scheme 
to treat such complex unsteady shock reflection problems as well as to provide detailed 
information on the flow physics of these flowfields. Experimental studies can provide useful 
information about a flow, but they can not supply with very detailed information. Therefore, 
the compliment of numerical studies is deemed important. 

Numerical Scheme 

The full Navier- Stokes equations in conservation law form are used in the present study: 


d t u + d x E + d v F = Re- l {d x E v + d y F v ) 
where U is a vector containing the conservation variables: 


U = 


fM 
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The vectors E and F are the inviscid flux vectors: 
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are the viscous flux vectors: 
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where r xx = (A -f 2 p)u x + Xv y) r xy — p{^y + v x), an d r yy — (^ + 2 p)v y -f Xu x . 

Define F ( = S^F + S^G, etc., where S ix and S iy are components of the surface vector 
of the numerical cells, then the Navier- Stokes equations can be discretized using a finite 
volume scheme as 


. +(F e - - (F< - (0.2) 

+(K) — — (^ij “ (0-3) 

= 0, (0.4) 

where v is the volume of the numerical cell. Thus the equation has been discretized in a 
general curvilinear coordinate system. The differences between various numerical schemes 
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are in the ways the fluxes on the numerical cell surfaces are constructed. It is well known that 
central difference causes odd and even decoupling and the dispersive error thus generated 
produces oscillations the the solution. While adding artificial viscosity can alleviate this 
problem, it often causes too much undesired numerical dissipation and hence the decadence 
of the solution. The TVD schemes developed in recent years are devised to overcome this 
difficulty by controlling the amount of artificial dissipation in the numerical formulation. 

In recent years, upwind schemes are generally preferred in solving hyperbolic equations 
because of the justifications found in theory of characteristics. In order to apply upwind 
difference or upwind interpolation, we need to first decompose an inviscid flux into the sum 
of one flux with only positive eigenvalues and one flux with only negative eigenvalues. In 
the present work, Roe’s flux- difference splitting is employed to decompose the fluxes. 

Roe’s flux-difference splitting is constructed by relating the flux difference at two ar- 
bitrary states, L and R, with the difference of variables U such that 


A F = AAU, 

A = A(U), A(U) = 

U = U(U L ,U R ). 

The object is to find an average state U so that the first equation is satisfied exactly for 
all admissible pairs ( Ul,Ur ). Since the eigenvalues of the Jacobian matrix A are real, the 
splitting of these eigenvalues according to their signs leads to the splitting of A: 


A = A + + A~. 
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The flux difference is split accordingly, 


A F = A F + + A F~ = A + AU + A~AU. 


Although the Roe splitting in Cartesian coordinates can be easily constructed, it is less 
straightforward in general coordinates within the framework of finite-volume discretization, 
where the geometric metrics are involved. A detailed derivation of the Roe splitting in 
three-dimensional general curvilinear coordinates was given by the present authors in Ref. 
1. 

Once the positive and negative flux differences are given, upwind interpolations for the 
positive and negative fluxes of various orders of accuracy can be constructed as follows. Let 



where the '+’ and components respectively correspond to nonnegative and nonpositive 
characteristic speeds. Thus upwind polynomial interpolation of fluxes using nodal (cell 


center) values gives the following formulas for the ‘positive’ flux: 

F+ +l = F? + dat{^-{l-K)A-F? + ai(l+K)A + F? (0.5) 

+a(3i(A + F t + - -A-l^ + )} = Ft + 6+, (0.6) 

<Xi J 

and for the ‘negative’ flux: 

F.“l = FCn ~ ^af+ijai+lt 1 “ (°- 7 ) 

+<r/3 i+1 (A + F:; 1 - — A-FS,)} = - «.-+,• (0-8) 

OLi+l J 
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where 


* = 0 

first-order upwind | 

9 = 1 

<7 — 0 

K= 1 

second-order central 

*= -1 

second-order full-upwind 

* = 0 

FVomm 

<7*1 

* = 0 

third-order biased 


and let l t — width of cell i in the ^-direction; 

4 + 4-i 

4 + 4+1 ’ 

4 

4 + 4+i ’ 

4 

(1 ± «)(/, + 4-i) + (1 T «)(4 + 4+i)' 

The first term in (2.3a) and (2.3b) is identified as the first-order upwind scheme, while 
the remaining higher-order terms are called the anti-diffusive terms [23]. Combining eqns 
(2.3a) and (2.3b), we obtain 

+ F m ) - i - F-) + («+ - S- +1 ). 

The first term is a simple-average representation of the interface flux [13], but neglects the 
nonuniformity of cells. The first two terms together constitute the first-order upwind scheme. 

It is a well known fact that first order upwind scheme, although very stable, is too 
diffusive. On the other hand, the dispersive errors contained in the higher order schemes 
cause oscillations in the solution. A compromise can be reached by using a limiter to limit 
the anti-diffusinve terms in eqs. (3) and (4) such that when oscillations occurs, the limiter 
would turn off the anti-diffusive terms; while for smooth solution, the formulas return to 
their higher order forms. 


A 
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The construction of a limited flux is nonlinear and requires modification in the high- 
order interpolation. This amounts to performing the following modifications in eqs. (3) and 
(4): 

F+ = Fi + +9at\±(l-K)<p{rt)&-F? 

2 l * 

+ ai (l + /c)*>(£)A+J* 

t 

+c0M7f'l^ +F * - i^)A-^ + ]}. 

= F,~+\ - «“.+i{°i+i(l - ' 1 )+( r ."+i) A * F .+i 

+^t( 1 + + 


The ‘limiters’ <p + and <p~ are functions of some appropriate ratio of neighboring flux 
differences. For a system of equations, there is no unique way of defining these functions. We 
have found the following set of definitions quite satisfactory over a wide range of problems 

[4,6,18]: 


+ _ <A + 

r i — <A-fl+,A-Ft> 
- _ 

r i ~ <a+f~,a+f-> 


The Roe’s ‘superbee limiter is used in the present study: 

tp(r) = mai{0,mtn(2r, l),mira(r,2)}. 


Time Evolution. The time integration integration scheme assembles fluxes entering 
the cell interfaces and updates the conservative variables at the cell center. The requirement 
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of a monotonic evolution places restrictions on the recombination stage and thus fixes the 
limiter function. As stated earlier, time accuracy is required in this evolution stage for 

transient problems. The present method originates from the Taylor series expansion in time, 
as was done in the Lax-Wendroff scheme. A two-step scheme with second order time accuracy 
is given here: 


( 1 ) : U* = U n + At2£, 

( 2 ) : = + 

U n+ 1 = \{U n + U**). 

The residual vanishes as the solution approaches a steady state. This scheme is similar in 
form to various types of predictor-corrector schemes. It is also restricted by the usual CFL 
stability condition. However, unlike the Lax-Wendroff two-step method, both predictor and 
corrector steps are identical, with no need of defining a midpoint for the corrector step. 
This reduces the complexity of evaluating the transport terms in the general curvilinear 
coordinates. The method recovers the one-step Lax-Wendroff scheme even for nonlinear flux 
functions, thus preserving second-order time accuracy. 


Results and Discussion 

Recently the problem of a blast wave passing an obstacle has attracted the attention 
of both experimentalists and numerical analysists. When a shock wave passes a circular 
cylinder, a serious of shock-solid boundary interaction occurs. The interaction produces an 
interesting shock reflection pattern. The unsteady nature and the intricate flow pattern made 
the problem a challenging one on which to test the ability of a shock capturing unsteady 
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flow solver. 


Figure 1 shows an experimental schlieren photograph given by Ref. 2. The experiment 
shows that when a planar shock wave first impinge on the cylinder, the shock has a regular 
reflection on the surface. As the reflection point moves towards downstream, it soon become 
a Mach reflection and the shock split into three, joined by a triple point. A slip surface 
(density discontinuity) forms at the triple point. After the moving shock has passed the 
cylinder, there remain a bow shock in front of the cylinder, and vortices are generated 
behind the cylinder. 

In the present work, a traveling shock of incident Mach number 2.81 impinging on 
a circular cylinder is numerically simulated with the recently developed upwind unsteady 
flow solver. The air surrounding the cylinder is originally at rest. The shock is released 
immediately in front of the cylinder at time equal to zero. The grid used is a polar coordinate 
system, with 100 points in the radial direction and 360 points in the circumferencial direction. 
The grid has a stretching in the radial direction. 

Figure 2 shows the flowfield properties as the shock just hit the cylinder. The density 
and pressure contours clearly show a regular shock reflection. In front of the cylinder, there 
exist a region of very high pressure; this can be observed from the pressure distribution on 
the cylinder surface, as shown in Figure 2c. This initial period produces the maximum load 
on the cylinder. It is also worth noting that, as shown in the surface pressure distribution, 
the reflection point (including two shocks) is captured within five grid points by the present 
numerical scheme. 

As the incident shock moved downstream, the shock reflection became a so called 
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Mach reflection. In another word, the intersection between the incident shock and the bow 
shock moved away from the solid surface, forming a triple point. A slip surface (density 
discontinuity) can be observed at this triple point from the density contours shown in Figure 
3a as wiggles in the contours. Though the slip surface does not have the sharp feature as 
observed in the schlieren photo. 

As shown in Figure 4, when the two branches of the incident shock collided at the rear 
of the cylinder, the crossing of the shocks created a local high pressure. A pair of week shocks 
can also be observed at about 25 degree and 335 degree; these week shocks are presumably 
caused by the flow over expansion behind the cylinder. 

When the incident shock moved further away from the cylinder, an interesting flow 
pattern is generated. The shocks seem branched out before they reach the solid surface. A 
pair of vortices appeared behind the cylinder. This flow pattern can be observed from the 
contours given in Figure 5. 


Concluding Remarks 

The newly developed upwind unsteady Navier-Stokes solver was used to study the 
problem of a blast wave interacting with a circular cylinder. The results show that the current 
solver can capture moving shocks quite well; it keeps the sharp feature of the incident shock 
as well as the reflected shocks. The slip surface or density can also be captured, although 
the calculated results do not show the sharp feature as seen on a schlieren photo. In general, 
one can conclude that moving shock problems can be numerically simulated with fairly good 
accuracy using an upwind scheme. 
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Figure 1. Schlieren photographs of shock wave reflection by a circular 
cylinder given by Ref. 2. 

(Reprinted with the permission of Cambridge University Press, 
Journal of Fluid Mechanics , Vol . 10, 1961, pp. 1-16, plates 
2 and 3, "Diffraction of Strong Shocks by Cones, Cylinders 
and Spheres" by A. E. Bryson and R. Gross) 
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Figure 2. Numerical solution of shock wave reflection: simple reflection; 
(a) density contours, (b) pressure contours, (c) Mach number contours, (d) 
wall pressure distribution. 
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Figure 3. Numerical solution of shock wave reflection: Mach reflection; 
(a) density contours, (b) pressure contours, (c) Mach number contours, (d) 
wall pressure distribution. 






Figure 4. Numerical solution of shock 
wave reflection: incident shock has just 
passed the cylinder; (a) density contours, 
(b) pressure contours, (c) Mach number 
contours, (d) wall pressure distribution. 
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